When not to avoid inbreeding: a gene’s eye view perspective
Author
Thomas Keaney, Arvid Agren and Hanna Kokko
Load packages
Code
# for tidy style coding and plottinglibrary(tidyverse) library(vroom) # to read lots of csv files at once# more table optionslibrary(pander) # for tableslibrary(kableExtra) # for scrolling tableslibrary(data.table) # for efficient handling of large dataframes# making ggplot more powerfullibrary(MetBrewer) # for colour palettes based upon artwork housed at the METlibrary(MoMAColors) # for colour palettes based upon artwork housed at MoMAlibrary(wesanderson) # for colour palettes based on wes anderson movieslibrary(PNWColors) # for colour palettes library(tidybayes) # for plotting distributionslibrary(stickylabeller) # labelling facets with strings in ggplotlibrary(geomtextpath) # for curved plot annotationslibrary(ggtext) # for markdown syntax in plot labelslibrary(patchwork) # for patching plots togetherlibrary(ggnewscale) # to reset scales in plots, allowing multiple fill arguments in ggplot# for computation speed checkslibrary(profvis) # breakdown of complex functionslibrary(bench) # individual functions
The model
The seminal equation
To model the inclusive fitness gained from an inbred mating, three components that contribute to fitness are required:
The number of offspring produced directly: \(n\)
The reduction in fitness (number of offspring) due to inbreeding: \(-\delta n\)
The indirect fitness gain (number of offspring) due to inbreeding: \(rn\), where \(r\) is the relatedness coefficient
Put together, the inclusive fitness from a single inbred mating is:
\[(1 + r)(1 - \delta)n\]
while fitness from a single outbred mating is simply \(n\).
When \((1 + r)(1 - \delta)n \gt n\) selection should favour a preference for inbreeding.
Solving the inequality for \(\delta\):
\[\delta \lt \frac{r}{1 + r}\] which for varying values of \(r\) looks like this:
Code
inbreeding_maximum_function <-function(r){ r / (1+ r)}parameters <-expand_grid(r =seq(from =0, to =1, by =0.05),delta =seq(from =0, to =1, by =0.05))r <- parameters %>%distinct(r)inbreeding_equilibria <-map_dfr(r, inbreeding_maximum_function) %>%rename(depression_threshold = r) %>%bind_cols(r)inbreeding_equilibria %>%ggplot(aes(x = r, y = depression_threshold)) +geom_line(linewidth =0.8) +coord_cartesian(ylim =c(0, 1)) +labs(x ='_r_, the relatedness coefficient',y =~delta~'(inbreeding depression)') +scale_x_continuous(expand =c(0, 0.009)) +scale_y_continuous(expand =c(0, 0)) +theme_bw() +theme(text =element_text(size =14),axis.title.x =element_markdown())
Code
# (prop fitness lost)\n that can be tolerated"
The parameter space above the curve shows where inbreeding avoidance should evolve, while the parameter space below the curve shows where inbreeding preference should evolve.
Accounting for sex differences in genetic architecture
As stated above, inclusive fitness in the absence of inbreeding depression is \((1 + r)n\). Here \(r\) represents the correlation between genotypes carried by interacting females and males, under the implicit assumption that loci appear equally in both sexes. However, given that there is sexual dimorphism in genetic architecture for many taxa, \(r\) does not sufficiently represent the correlation between genotypes for all loci.
To delineate differences in the effect of \(r\) for different regions of the genome, we multiply \(r\) with a new variable \(a\), the probability that a locus present in one sex is also present in the gametes produced by the other. Unlike \(r\) which is relative to the population mean relatedness, \(a\) is expressed as an absolute value ranging from 0 to 1.
The indirect component of fitness accrued by from an inbred mating becomes
\[ran\]
and inclusive fitness from an inbreeding event becomes
\[(1 + ra)(1 - \delta )n\]
Taking an allele found at a diploid autosomal locus as an example, all of the gametes produced by a relative possess this locus, where they could potentially carry alleles identical by descent. In this case \(a = 1\) and the indirect component of inclusive fitness is dictated solely by \(r\). The results for this autosomal scenario are presented in Parker (1979), Kokko and Ots (2006) and others who have explored this topic. In contrast, an inbreeding preference allele present at a locus on a Y or W chromosome has no opportunity to propagate any alleles identical by descent through inbreeding, as these chromosomes are not carried by the gametes of the opposite sex mating partner. In this case \(a = 0\). However, as inbreeding depression is a result of homozygosity for deleterious recessive alleles throughout the genome, the costs of inbreeding depression are born by all alleles carried by the individual. Conflict over the expression of inbreeding preference between alleles present on autosomes and those present on hemizygous sex chromosomes is immediately clear.
X- or Z-linked loci present an interesting intermediate case, with sex-specific values for \(a\). When the inbreeding locus is carried by the sex with homozygous sex chromosomes, \(a\) is half that of autosomal loci, whereas it does not depart from the autosomal case when the locus is found within the hemizygous sex. Using loci on the X as an example, those present in a XX female are only found in ~50% of a interacting males gametes, as the remaining 50% carry Y chromosomes (assuming an even primary sex ratio). When an X-linked locus is found in a male, an interacting female’s gametes all carry X chromosomes and \(a = 1\).
The X/Z situation is made additionally complex because there is an element of frequency dependence to the kin selected benefits. When an inbreeding allele on an autosome is rare, then the chance of a relative carrying two copies is low, whereas when the allele is common, this chance is much higher. Rarity therefore leads to similar fitness outcomes for autosomal and X/Z linked alleles (when present in the hemizygous sex), while commonality of the allele likely roughly equates to the conflicting situation outlined in the above paragraph. However, frequency dependence might not as relevant as I initially expected, because the inbreeding allele becomes very common quickly within families, even whilst rare across the population (first proposed in Fisher, 1930).
Table 1. Values of the parameter \(a\) for different regions of the genome. \(a\) is the one-way probability that a locus carried by one individual is found within the gametes of an opposite sex individual. Note that cytoplasmic chromosomes are assumed to have exclusive maternal inheritance.
Code
x <-c(1, # autosomes, X chromosome males or Z chromosome females, haplodiploid both sexes when producing females 0, # Y or W chromosome0.5# X chromosome females or Z chromosome males )tibble(`Prob. that opposite sex gametes carry focal locus`=c(1, 0.5, 0),`Relevant cases`=c("Autosomes in either sex, X chromosomes in males, Z chromosomes in females, chromosomes in haplodiploids of either sex when reproducing sexually, cytoplasmic chromosomes in males","X chromosomes in females, Z chromosomes in males","Y chromosomes in males, W chromosomes in females, cytoplasmic chromosomes in females")) %>%pander(split.cell =20, split.table =Inf)
Prob. that opposite sex gametes carry focal locus
Relevant cases
1
Autosomes in either sex, X chromosomes in males, Z chromosomes in females, chromosomes in haplodiploids of either sex when reproducing sexually, cytoplasmic chromosomes in males
0.5
X chromosomes in females, Z chromosomes in males
0
Y chromosomes in males, W chromosomes in females, cytoplasmic chromosomes in females
Once again we can find the condition where breeding with a relative returns greater fitness than an inbreeding avoidance strategy, this time accounting for genetic architecture
\[\delta \lt \frac{ra}{1 + ra}\] Ignoring frequency dependence for now, we can plot the new slopes produced by varying \(r\) and \(a\):
Differences between the sexes beyond genetic architecture
Parker’s seminal equations:
In his 1979 book chapter, Parker considered the inclusive fitness results of breeding with a relative and identified that females and males should have different tolerances for inbreeding depression. The key departure from the unlimited polygyny case presented in the above equations is that a cost to future reproductive success is included for males, i.e. due to finite sperm production, parental care, or harmful mating behaviour such as sexual cannibalism.
For males, Parker found that selection would favour inbreeding with a sister (full-sib) who could otherwise outcross when:
\[n(1 - \delta) + rn(1- \delta) - cn \gt rn\]
the first term is the direct number of alleles propagated, the second term is the indirect number of alleles propagated (note that this is weighted by relatedness), the third term is the direct number of alleles that were not directly propagated by the male through outcrossing, and the final opposing term is the number of alleles that would’ve been transmitted had his sister outcrossed (and he forgone mating).
\(c\) is the cost of the present mating, relative to what is lost for a female. This can be considered a ratio of parental investment. When \(c = 1\) parental investment in the current bout of reproduction is even between the sexes. Alternatively, if males only contribute cheaply produced sperm to an incestuous mating, the cost of mating is likely very small relative to females i.e. \(c ~ 0\).
We add the \(a\) variable to the equation and letting \(n = 1\), simplify to
\[(1-\delta) + ra(1-\delta) - c \gt ra\]
We can again find the condition where breeding with a relative returns greater fitness than an inbreeding avoidance strategy:
\[\delta_{male} = \frac{1 - c}{1 + ra}\]
Parker then modelled the condition for monandrous females to prefer incestuous matings when also presented with an outcrossing opportunity.
\[n(1 - \delta) + rn(1-\delta) - crn \gt n\]
which we can write as
\[(1-\delta) + ra(1-\delta) - rac \gt 1\]
the inbreeding depression threshold is
\[\delta_{female} = \frac{ra - rac}{1 + ra}\]
Note that when \(c = 0\), this is equivalent to the \(\delta\) threshold found in the single mating case.
Plot the relationship between \(r\) and \(\delta\) for several values of \(c\) and \(a\)
In the figure below, the dashed lines indicate the level of inbreeding depression that can be tolerated for a given value of \(r\). The plot is split into panels by \(c\), the cost of mating for males relative to females and \(a\), the probability that a locus present in one sex is also present in the gametes produced by the other.
Code
pal1 <-met.brewer("OKeeffe1", n=100, direction =-1)pal2 <-met.brewer("Hiroshige", n=50, direction =-1)#pal2 <- moma.colors("Avedon", n = 50, direction = 1)oranges <-c("#ffe6b7", "#ffd06f", "#f7aa58", "#ef8a47", "#e76254") blues <-c("#aadce0", "#72bcd5", "#528fad", "#376795", "#1e466e")# Lowest colour for all gradientsgradient_base <- oranges[1]gradient_max <- oranges[5]# Create a list of gradients for each colour 2 to 10 over five steps from # gradient_base grey (low) to colour (high)my_gradient <-colorRampPalette(c(gradient_base, gradient_max))(50)Female_plot <- analytical_results %>%filter(c ==0| c ==0.5) %>%ggplot(aes(x = r, y = D)) +geom_blank() +geom_raster(aes(fill = female_fitness_contrast)) +stat_contour(aes(z = female_fitness_contrast), colour ="black", binwidth =0.25,breaks =c(-1.25, -1, -0.75, -0.5, -0.25, 0.25, 0.5, 0.75, 1, 1.25)) +stat_contour(aes(z = female_fitness_contrast), colour ="black", breaks =0,linetype =2) +scale_fill_gradientn(colours = pal1, breaks =c(-2, -1, 0, 1, 2), limits =c(-2, 2)) +facet_wrap(c ~ a, scales ="free", nrow =2, strip.position =c("top"),labeller =label_glue('Prob. related gamete carries inbreeding locus = {`a`}\nMale mating cost = {`c`}')) +labs(x ='_r_, the relatedness coefficient',y =~delta~'(inbreeding depression)',fill ="Inbreeding fitness",title ="A. Alleles present in females") +scale_x_continuous(expand =c(0, 0)) +scale_y_continuous(expand =c(0, 0)) +# labels = c(0, 25, 50, 75, 90)) +theme(panel.border =element_rect(fill =NA, colour ="black", linewidth = .8),strip.background =element_rect(colour ="black", fill ="Aliceblue", linewidth = .8),strip.text =element_text(size =7),axis.title.x =element_markdown(),plot.title =element_text(hjust =0.5))Male_plot <- analytical_results %>%filter(c ==0| c ==0.5) %>%ggplot(aes(x = r, y = D)) +geom_blank() +geom_raster(aes(fill = male_fitness_contrast)) +stat_contour(aes(z = male_fitness_contrast), colour ="black", binwidth =0.25,breaks =c(-1.25, -1, -0.75, -0.5, -0.25, 0.25, 0.5, 0.75, 1, 1.25)) +stat_contour(aes(z = male_fitness_contrast), colour ="black", breaks =0,linetype =2) +scale_fill_gradientn(colours = pal1, breaks =c(-2, -1, 0, 1, 2), limits =c(-2, 2)) +facet_wrap(c ~ a, scales ="free", nrow =2, strip.position =c("top"),labeller =label_glue('Prob. related gamete carries inbreeding locus = {`a`}\nMale mating cost = {`c`}')) +labs(x ='_r_, the relatedness coefficient',y =~delta~'(inbreeding depression)',fill ="Inbreeding fitness",title ="B. Alleles present in males") +scale_x_continuous(expand =c(0, 0)) +scale_y_continuous(expand =c(0, 0)) +# labels = c(0, 25, 50, 75, 90)) +theme(panel.border =element_rect(fill =NA, colour ="black", linewidth = .8),strip.background =element_rect(colour ="black", fill ="Aliceblue", linewidth = .8),strip.text =element_text(size =7),axis.title.x =element_markdown(),plot.title =element_text(hjust =0.5))Female_plot / Male_plot +plot_layout(guides ="collect")
When is there evolutionary conflict over inbreeding?
To quantify conflict in its various forms we use a version of the I index presented in Innocenti and Morrow (2010).
Code
Intragenomic_conflict_females <- analytical_results %>%select(1:4, female_fitness_contrast) %>%pivot_wider(names_from = a, values_from = female_fitness_contrast) %>%#mutate(`Autosome - X` = `1` - `0.5`,# `Autosome - Z` = 0,# `Autosome - W & Z - W` = `1` - `0`) %>% mutate(`Autosome vs X`=`1`*`0.5`/sqrt(((`1`)^2+ (`0.5`)^2)/2),`Autosome vs Z`=`1`*`1`/sqrt(((`1`)^2+ (`1`)^2)/2),`Autosome & Z vs W`=`1`*`0`/sqrt(((`1`)^2+ (`0`)^2)/2)) %>%pivot_longer(cols =7:9, names_to ="contrast", values_to ="Evolutionary_conflict") %>%mutate(relationship =case_when( contrast =="Autosome vs X"&`1`>0&`0.5`>0~"Inbreeding favoured in both contexts", contrast =="Autosome vs X"&`1`<0&`0.5`<0~"Inbreeding deleterious in both contexts", contrast =="Autosome vs Z"&`1`>0~"Inbreeding favoured in both contexts", contrast =="Autosome vs Z"&`1`<0~"Inbreeding deleterious in both contexts", contrast =="Autosome & Z vs W"&`1`>0&`0`>0~"Inbreeding favoured in both contexts", contrast =="Autosome & Z vs W"&`1`<0&`0`<0~"Inbreeding deleterious in both contexts",`0`==0&`0.5`==0&`1`==0~"Inbreeding deleterious in both contexts", .default ="Intragenomic conflict")) %>%mutate(contrast =fct_relevel(contrast, "Autosome vs X", "Autosome vs Z","Autosome & Z vs W")) %>%filter(c ==0) # remove if we want more c valuesIntragenomic_conflict_males <- analytical_results %>%select(1:4, male_fitness_contrast) %>%pivot_wider(names_from = a, values_from = male_fitness_contrast) %>%#mutate(`Autosome - X` = 0,# `Autosome - Z` = `1` - `0.5`,# `Autosome - Y & X - Y` = `1` - `0`) %>%mutate(`X vs Autosome`=`1`*`1`/sqrt(((`1`)^2+ (`1`)^2)/2),`Z vs Autosome`=`1`*`0.5`/sqrt(((`1`)^2+ (`0.5`)^2)/2),`Y vs Autosome & X`=`1`*`0`/sqrt(((`1`)^2+ (`0`)^2)/2)) %>%pivot_longer(cols =7:9, names_to ="contrast", values_to ="Evolutionary_conflict") %>%mutate(relationship =case_when( contrast =="X vs Autosome"&`1`>0~"Inbreeding favoured in both contexts", contrast =="X vs Autosome"&`1`<0~"Inbreeding deleterious in both contexts", contrast =="Z vs Autosome"&`1`>0&`0.5`>0~"Inbreeding favoured in both contexts", contrast =="Z vs Autosome"&`1`<0&`0.5`<0~"Inbreeding deleterious in both contexts", contrast =="Y vs Autosome & X"&`1`>0&`0`>0~"Inbreeding favoured in both contexts", contrast =="Y vs Autosome & X"&`1`<0&`0`<0~"Inbreeding deleterious in both contexts",`0`==0&`0.5`==0&`1`==0~"Inbreeding deleterious in both contexts", .default ="Intragenomic conflict")) %>%mutate(contrast =fct_relevel(contrast, "X vs Autosome", "Z vs Autosome","Y vs Autosome & X")) %>%filter(c ==0) # remove if we want more c valuesmake_genomic_conflict_plot <-function(data, enter_title, colour_pal){ data %>%ggplot(aes(x = r, y = D)) +geom_blank() +geom_tile(data = data %>%filter(relationship =="Intragenomic conflict"),aes(fill = Evolutionary_conflict*-1)) +#geom_tile(data = data,# aes(fill = Evolutionary_conflict)) + scale_fill_gradientn(colours = colour_pal, limits =c(0, 0.6), #na.value = "white",breaks =c(0, 0.6),labels =c("Weaker conflict", "Stronger conflict")) +labs(fill ="Evolutionary conflict") +new_scale_fill() +geom_tile(data = data %>%filter(relationship !="Intragenomic conflict"),aes(fill = relationship), alpha =0.75) +scale_fill_manual(values =c("#fbe6c5", "#d2fbd4"), labels =c("Inbreeding deleterious\nin both contexts", "Inbreeding favoured\nin both contexts")) +stat_contour(aes(z = Evolutionary_conflict*-1), colour ="black", binwidth =25,breaks =c(0, 0.2, 0.4, 0.6)) +facet_wrap(~contrast, nrow =3,scales ="free", strip.position =c("top"),labeller =label_glue('{`contrast`}')) +labs(x ='_r_ (relatedness coefficient)',y =~delta~'(inbreeding depression)',fill ="Evolutionary concordance",title = enter_title) +scale_x_continuous(expand =c(0, 0)) +scale_y_continuous(expand =c(0, 0)) +# labels = c(0, 25, 50, 75, 90)) +theme(panel.border =element_rect(fill =NA, colour ="black", size = .8),panel.grid.minor =element_blank(),strip.background =element_rect(colour ="black", fill ="Aliceblue", linewidth = .8),axis.title.x =element_markdown(),plot.title =element_text(hjust =0.5, size =7)) }icf <-make_genomic_conflict_plot(Intragenomic_conflict_females, "A. Female inter-chromosomal conflict", oranges)icm <-make_genomic_conflict_plot(Intragenomic_conflict_males, "B. Male inter-chromosomal conflict", oranges)
Make column 3
Code
autosomal_data <- analytical_results %>%filter(a ==1, c !=1, c !=0.25) %>%#mutate(Evolutionary_conflict = male_fitness_contrast - female_fitness_contrast) %>%mutate(Evolutionary_conflict = female_fitness_contrast * male_fitness_contrast /sqrt(((female_fitness_contrast)^2+ (male_fitness_contrast)^2)/2)) %>%mutate(relationship =case_when(female_fitness_contrast <0& male_fitness_contrast >0~"Sexual conflict", female_fitness_contrast <0& male_fitness_contrast <0~"Inbreeding deleterious in both contexts", female_fitness_contrast >0& male_fitness_contrast >0~"Inbreeding favoured in both contexts"),Location ="Autosome") X_data <- analytical_results %>%filter(a ==0.5, c !=1, c !=0.25) %>%select(1:4, contains("female")) %>%rename(a_female = a) %>%# this step makes the join work as intendedleft_join( analytical_results %>%filter(a ==1, c !=1, c !=0.25) %>%select(1:4, starts_with("male")) %>%rename(a_male = a) # this step makes the join work as intended ) %>%#mutate(Evolutionary_conflict = male_fitness_contrast - female_fitness_contrast) %>% mutate(Evolutionary_conflict = female_fitness_contrast * male_fitness_contrast /sqrt(((female_fitness_contrast)^2+ (male_fitness_contrast)^2)/2)) %>%mutate(relationship =case_when(female_fitness_contrast <0& male_fitness_contrast >0~"Sexual conflict", female_fitness_contrast <0& male_fitness_contrast <0~"Inbreeding deleterious in both contexts", female_fitness_contrast >0& male_fitness_contrast >0~"Inbreeding favoured in both contexts"),Location ="X")Z_data <- analytical_results %>%filter(a ==1, c !=1, c !=0.25) %>%select(1:4, contains("female")) %>%rename(a_female = a) %>%left_join( analytical_results %>%filter(a ==0.5, c !=1, c !=0.25) %>%select(1:4, starts_with("male")) %>%rename(a_male = a) ) %>%#mutate(Evolutionary_conflict = male_fitness_contrast - female_fitness_contrast) %>%mutate(Evolutionary_conflict = female_fitness_contrast * male_fitness_contrast /sqrt(((female_fitness_contrast)^2+ (male_fitness_contrast)^2)/2)) %>%mutate(relationship =case_when(female_fitness_contrast <0& male_fitness_contrast >0~"Sexual conflict", female_fitness_contrast <0& male_fitness_contrast <0~"Inbreeding deleterious in both contexts", female_fitness_contrast >0& male_fitness_contrast >0~"Inbreeding favoured in both contexts"),Location ="Z")plotting_data <-bind_rows(autosomal_data, X_data, Z_data) %>%filter(c ==0)Sexual_conflict_plot <- plotting_data %>%ggplot(aes(x = r, y = D)) +geom_blank() +geom_tile(data = plotting_data %>%filter(relationship =="Sexual conflict"),aes(fill = Evolutionary_conflict*-1)) +#scale_fill_gradientn(colours = pal2, limits = c(-1.2, 1.1), #na.value = "white",# labels = c("Strong conflict, female (+)", -0.5, # "No conflict", 0.5, "Strong conflict, male (+)")) +scale_fill_gradientn(colours = oranges, limits =c(0, 0.6), #na.value = "white",breaks =c(0, 0.6),labels =c("Weaker conflict", "Stronger conflict")) +labs(fill ="Evolutionary conflict") +new_scale_fill() +geom_tile(data = plotting_data %>%filter(relationship !="Sexual conflict"),aes(fill = relationship), alpha =0.75) +scale_fill_manual(values =c("#fbe6c5", "#d2fbd4"), labels =c("Inbreeding deleterious\nin both contexts", "Inbreeding favoured\nin both contexts")) +stat_contour(aes(z = Evolutionary_conflict*-1), colour ="black", binwidth =25,breaks =c(0, 0.2, 0.4, 0.6)) +facet_wrap(~Location, nrow =3,scales ="free", strip.position =c("top"),labeller =label_glue('{`Location`}')) +labs(x ='_r_ (relatedness coefficient)',y =~delta~'(inbreeding depression)',fill ="Evolutionary concordance",title ="C. Intra-chromosomal sexual conflict") +scale_x_continuous(expand =c(0, 0)) +scale_y_continuous(expand =c(0, 0)) +# labels = c(0, 25, 50, 75, 90)) +theme(panel.border =element_rect(fill =NA, colour ="black", size = .8),strip.background =element_rect(colour ="black", fill ="Aliceblue", linewidth = .8),axis.title.x =element_markdown(),plot.title =element_text(hjust =0.5, size =7))
Panels in columns A and B show the parameter space where inbreeding increases the propagation of one type of chromosome, but has an opposite, deleterious affect on the propagation of a second chromosome class - so called inter-chromosomal conflict.
Column A shows regions of inter-chromosomal conflict within females, where chromosome classes with a high probability of encountering gametes carrying the inbreeding locus have high fitness in the conflict zone. In females, autosomal and Z chromosomes are aligned in the fitness outcomes of inbreeding; these cases reflect the implicit assumptions of previous models exploring the inclusive fitness effects of inbreeding. In some cases, the fitness interests of chromosome classes align e.g. autosomes and Z chromosomes (this is only true for females).
The second column shows inter-chromosomal in males. Importantly, the direction of selection for inbreeding is reversed in comparison to the female case in the zones of conflict. In males, chromosome classes with high probabilities of inbreeding loci being found in female gametes now have low fitness in the conflict space.
Panels in column C show where loci are expected to be under intra-chromosomal sexual conflict over inbreeding preference. In regions of sexual conflict, inbreeding preference is always favoured in males, but has negative fitness consequences if expressed by females (assuming that males invest less into mating than females). Note that intra-chromosomal conflict encompasses both the intra- and inter-locus forms of sexual conflict.
Make Figure 2
We can also plot the cost to female direct reproductive output,
Figure 2. direct female reproductive output at the threshold where chromosome-specific selection should no longer favour inbreeding. Reproductive output is expressed relative to females that exclusively outcross.
The simulation
Parker’s \(c\) value offers a simple, intuitive way to model the cost of mating for males relative to females. However, other methods better capture the dynamics of real populations, where both sexes also run the risk of going unmated. To incorporate both costs of mating and matelessness, we simulate the invasion of an allele that encodes a preference for inbreeding, for loci on various chromosomes that are expressed in either females or males.
The simulation progresses through time by jumping from event to event. Events are calculated via a gillespie-like algorithm, which involves both memory-less processes and others that are maintained through time. Events can trigger state changes for the individuals in the population, leading to death, mating, inheritance of a breeding site or offspring production. The simulation ends when the inbreeding allele is purged from the population, it reaches a high frequency in the population, the population goes extinct, or the time limit elapses.
Following Ekrem and Kokko (2023), we find the lifespans of N diploid individuals in a sexually reproducing population to initialise the simulation. Mortality events are drawn from an exponential distribution with rate \(\lambda = 1/\mu\), where \(\mu\) is the mean time till death. The probability of mortality is therefore constant across life for each sex.
Lets have a look at the mean lifespan drawn from the exponential distribution parameterised with different mortality rates:
Code
#exponential_draws <-tibble(`Mean lifespan`=seq(from =0.001, to =1.001, by =0.1)) %>%mutate(`Corresponding mortality rate`=1/`Mean lifespan`) %>%arrange(-`Mean lifespan`) %>%pander()
Mean lifespan
Corresponding mortality rate
1.001
0.999
0.901
1.11
0.801
1.248
0.701
1.427
0.601
1.664
0.501
1.996
0.401
2.494
0.301
3.322
0.201
4.975
0.101
9.901
0.001
1000
We use these values to inform our parameterisation of inbreeding depression (\(\delta\)). Individuals that don’t inbreed produce offspring with a mean lifespan of 1 time unit (\(\lambda = 1\)). For inbred offspring, we model inbreeding depression by increasing the risk of mortality (\(\lambda > 1\)), lowering the mean lifespan of inbred individuals relative to their outbred counterparts.
Each individual possesses an inbreeding locus with two possible alleles. The A_O allele encodes inbreeding avoidance, whereas the A_I allele encodes inbreeding preference. The number_mutants argument allows us to set the number of mutations that the population is originally seeded with.
Female-male encounters are found using the mass action law (also called fertilisation kinetics in this biological context), with male search efficiency \(v\). To get an idea of how \(v\) influences mating encounters, we plot the population wide male-female encounter timestamps. The dashed line shows the mean lifespan of outbred individuals in our simulations. Note that these are just encounters; whether they lead to mating depends upon past mating experiences of both females and males.
Figure. SX: the proportion of the female population encountered by a single male across his lifetime, split by how efficiently males search for females. If there are 100 females in the population, a value of 0.1 indicates the male encounters 10 females on average. The decline in potential mating encounters as lifespan decreases indicates the severity of the inbreeding depression cost for males.
The cost of inbreeding depression manifests differently for females, where the ability to reproduce depends on mating with a male and securing a breeding site. If the population is at carrying capacity, the latter criterion is a function of age, as breeding sites are allocated randomly to non-breeding females once a breeder dies. Therefore, the longer a female lives, the more chances she has of becoming a breeder. If the population is below carrying capacity, the cost of inbreeding is less severe, as a new female can immediately secure a breeding site, leaving mating as the only restriction to offspring production.
Build a mating table
To simulate the invasion of alleles found on different chromosomes, we consider simple single-locus, two-allele dynamics. For each chromosomal inheritance system, we build a mating table, which tracks the production of new genotypes from each possible mother-father pair of genotypes in the population. Inheritance is mendelian, except for cytoplamsic chromosomes, which are exclusively maternally inherited.
The simulation tracks the invasion of an allele encoding a preference for inbreeding (denoted I), into a population initially dominated by an allele that encodes inbreeding avoidance (denoted O for outbreeding).
Code
make_mating_table <-function(gene_location){ make_offspring <-function(X, Y, type, zygote_freq, gene_location){tibble(Female_genotype = X,Male_genotype = Y, type, zygote_freq,locus_type = gene_location) }# Specify the possible offspring genotypes for all the potential crosses; we use these for the type argument in the make_offspring function# autosomal# II x II a_genotype_1 <-c("A_IA_I.Female", "A_IA_I.Male")# II x IO a_genotype_2 <-c("A_IA_I.Female", "A_IA_I.Male", "A_IA_O.Female", "A_IA_O.Male")# II x OO a_genotype_3 <-c("A_IA_O.Female", "A_IA_O.Male")# IO x IO a_genotype_4 <-c("A_IA_I.Female", "A_IA_I.Male", "A_IA_O.Female", "A_IA_O.Male", "A_OA_O.Female", "A_OA_O.Male")# IO x OO a_genotype_5 <-c("A_IA_O.Female", "A_IA_O.Male","A_OA_O.Female", "A_OA_O.Male")# OO x OO a_genotype_6 <-c("A_OA_O.Female", "A_OA_O.Male")# XY# II x IY_I xy_genotype_1 <-c("X_IX_I.Female", "X_IY_I.Male")# II x IY_O xy_genotype_2 <-c("X_IX_I.Female", "X_IY_O.Male")# II x OY_I xy_genotype_3 <-c("X_IX_O.Female", "X_IY_I.Male")# II x OY_O xy_genotype_4 <-c("X_IX_O.Female", "X_IY_O.Male")# IO x IY_I xy_genotype_5 <-c("X_IX_I.Female", "X_IY_I.Male","X_IX_O.Female", "X_OY_I.Male")# IO x IY_O xy_genotype_6 <-c("X_IX_I.Female", "X_IY_O.Male", "X_IX_O.Female", "X_OY_O.Male")# IO x OY_I xy_genotype_7 <-c("X_IX_O.Female", "X_IY_I.Male","X_OX_O.Female", "X_OY_I.Male")# IO x OY_O xy_genotype_8 <-c("X_IX_O.Female", "X_IY_O.Male","X_OX_O.Female", "X_OY_O.Male")# OO x IY_I xy_genotype_9 <-c("X_IX_O.Female", "X_OY_I.Male")# OO x IY_O xy_genotype_10 <-c("X_IX_O.Female", "X_OY_O.Male")# OO x OY_I xy_genotype_11 <-c("X_OX_O.Female", "X_OY_I.Male")# OO x OY_O xy_genotype_12 <-c("X_OX_O.Female", "X_OY_O.Male")# ZW# IW_I x II zw_genotype_1 <-c("Z_IZ_I.Male", "Z_IW_I.Female")# IW_I x IO zw_genotype_2 <-c("Z_IZ_I.Male", "Z_IZ_O.Male", "Z_IW_I.Female", "Z_OW_I.Male")# IW_I x OO zw_genotype_3 <-c("Z_IZ_O.Male", "Z_OW_I.Female")# IW_O x II zw_genotype_4 <-c("Z_IZ_I.Male", "Z_IW_O.Female")# IW_O x IO zw_genotype_5 <-c("Z_IZ_I.Male", "Z_IZ_O.Male","Z_IW_O.Female", "Z_OW_O.Female")# IW_O x OO zw_genotype_6 <-c("Z_IZ_O.Male", "Z_OW_O.Female")# OW_I X II zw_genotype_7 <-c("Z_IZ_O.Male", "Z_IW_I.Female")# OW_I x IO zw_genotype_8 <-c("Z_IZ_O.Male", "Z_OZ_O.Male","Z_IW_I.Female", "Z_OW_I.Female")# OW_I x OO zw_genotype_9 <-c("Z_OZ_O.Male", "Z_OW_I.Female")# OW_O X II zw_genotype_10 <-c("Z_IZ_O.Male", "Z_IW_O.Female")# OW_O x IO zw_genotype_11 <-c("Z_IZ_O.Male", "Z_OZ_O.Male","Z_IW_O.Female", "Z_OW_O.Female")# OW_O x OO zw_genotype_12 <-c("Z_OZ_O.Male", "Z_OW_O.Female")# cytoplasmic# I x I# I x O c_genotype_1 <-c("C_I.Female", "C_I.Male")# O x O# O x I c_genotype_2 <-c("C_O.Female", "C_O.Male")# Now calculate the zygote frequencies for each cross# autosomal# even frequency of two offspring genotypes freq_2 <-rep(0.5, 2)# even frequency between four offspring types freq_4 <-rep(0.25, 4)# when there are 6 offspring genotypes freq_6 <-c(0.125, 0.125,0.25, 0.25,0.125, 0.125)bind_rows(list(make_offspring("A_IA_I", "A_IA_I", a_genotype_1, freq_2, "autosomal"),make_offspring("A_IA_I", "A_IA_O", a_genotype_2, freq_4, "autosomal"),make_offspring("A_IA_I", "A_OA_O", a_genotype_3, freq_2, "autosomal"),make_offspring("A_IA_O", "A_IA_I", a_genotype_2, freq_4, "autosomal"),make_offspring("A_IA_O", "A_IA_O", a_genotype_4, freq_6, "autosomal"),make_offspring("A_IA_O", "A_OA_O", a_genotype_5, freq_4, "autosomal"),make_offspring("A_OA_O", "A_IA_I", a_genotype_3, freq_2, "autosomal"),make_offspring("A_OA_O", "A_IA_O", a_genotype_5, freq_4, "autosomal"),make_offspring("A_OA_O", "A_OA_O", a_genotype_6, freq_2, "autosomal"),make_offspring("X_IX_I", "X_IY_I", xy_genotype_1, freq_2, "XY"),make_offspring("X_IX_I", "X_IY_O", xy_genotype_2, freq_2, "XY"),make_offspring("X_IX_I", "X_OY_I", xy_genotype_3, freq_2, "XY"),make_offspring("X_IX_I", "X_OY_O", xy_genotype_4, freq_2, "XY"),make_offspring("X_IX_O", "X_IY_I", xy_genotype_5, freq_4, "XY"),make_offspring("X_IX_O", "X_IY_O", xy_genotype_6, freq_4, "XY"),make_offspring("X_IX_O", "X_OY_I", xy_genotype_7, freq_4, "XY"),make_offspring("X_IX_O", "X_OY_O", xy_genotype_8, freq_4, "XY"),make_offspring("X_OX_O", "X_IY_I", xy_genotype_9, freq_2, "XY"),make_offspring("X_OX_O", "X_IY_O", xy_genotype_10, freq_2, "XY"),make_offspring("X_OX_O", "X_OY_I", xy_genotype_11, freq_2, "XY"),make_offspring("X_OX_O", "X_OY_O", xy_genotype_12, freq_2, "XY"),make_offspring("Z_IW_I", "Z_IZ_I", zw_genotype_1, freq_2, "ZW"),make_offspring("Z_IW_I", "Z_IZ_O", zw_genotype_2, freq_4, "ZW"),make_offspring("Z_IW_I", "Z_OZ_O", zw_genotype_3, freq_2, "ZW"),make_offspring("Z_IW_O", "Z_IZ_I", zw_genotype_4, freq_2, "ZW"),make_offspring("Z_IW_O", "Z_IZ_O", zw_genotype_5, freq_4, "ZW"),make_offspring("Z_IW_O", "Z_OZ_O", zw_genotype_6, freq_2, "ZW"),make_offspring("Z_OW_I", "Z_IZ_I", zw_genotype_7, freq_2, "ZW"),make_offspring("Z_OW_I", "Z_IZ_O", zw_genotype_8, freq_4, "ZW"),make_offspring("Z_OW_I", "Z_OZ_O", zw_genotype_9, freq_2, "ZW"),make_offspring("Z_OW_O", "Z_IZ_I", zw_genotype_10, freq_2, "ZW"),make_offspring("Z_OW_O", "Z_IZ_O", zw_genotype_11, freq_4, "ZW"),make_offspring("Z_OW_O", "Z_OZ_O", zw_genotype_12, freq_2, "ZW"),make_offspring("C_I", "C_I", c_genotype_1, freq_2, "cytoplasmic"),make_offspring("C_I", "C_O", c_genotype_1, freq_2, "cytoplasmic"),make_offspring("C_O", "C_I", c_genotype_2, freq_2, "cytoplasmic"),make_offspring("C_O", "C_O", c_genotype_2, freq_2, "cytoplasmic") )) %>%filter(locus_type == gene_location)}
continuous_time_simulation <-function(row, parameters, inheritance_scheme){# package loading need to be done within the function to run on a cluster - local environment is not inherited by clusterlibrary(data.table) library(tidyverse) sample_vec <-function(x, ...) x[sample(length(x), ...)] # so we can sample from vectors with length 1 without this being interpreted as an integerprint(paste("Doing row", row)) # this shows which row in the parameter space is being modelled prop_i_table <-data.table(time =numeric(), proportion_I =numeric(),population_size =numeric()) # filled in as sim progresses keep_going <-TRUE# if the inbreeding allele fixes or goes extinct, this will change to false and the while loop will quit early Starting_pop_size <- parameters$Starting_pop_size[row] N <- parameters$N[row] # constant number_mutants <- parameters$number_mutants[row] # constant at 1 baseline_mean_lifespan <- parameters$baseline_mean_lifespan[row] # constant at 2 time_end <- parameters$time_end[row] # a cut-off point for each run sex_expressed <- parameters$sex_expressed[row] chromosome <- parameters$chromosome[row] heterozygous_genotype <- parameters$heterozygous_genotype[row] homozygous_genotype <- parameters$homozygous_genotype[row] hemizygous_genotype <- parameters$hemizygous_genotype[row]#C <- parameters$C[row] v <- parameters$v[row] refractory_period <- parameters$refractory_period[row] D <- parameters$D[row] dominance <- parameters$dominance[row] parameter_space_ID <- parameters$parameter_space_ID[row]# define the starting genotypes for each sex so the population table can be built Female_starting_genotype <- inheritance_scheme[.N]$Female_genotype Male_starting_genotype <- inheritance_scheme[.N]$Male_genotype# Set the number of breeding sites breeding_sites <-0.5*Starting_pop_size# Initialize the Individual_ID and Family_ID counters Individual_ID_counter <- Starting_pop_size Family_ID_counter <- Starting_pop_size/N # family size equals the no. offspring produced by a single female# the simulation tracks the population via a data.table# create the starting population - note that females are sex = 1 and males are sex = 0 population <-data.table(mortality_rate =1/baseline_mean_lifespan,Sex =rbinom(n = Starting_pop_size, 1, prob =0.5),birth_time =0,matings =0,reproduced =0,mated_with ="NA",inbred_mating =0,refractory_period_end =0,relative_encountered =0 )[, `:=` (Genotype =ifelse(Sex >0, Female_starting_genotype, Male_starting_genotype),breeding = Sex,Individual_ID = .I,Family_ID =rep(1:(.N/N), each = N, length.out = .N))]# seed population with the inbreeding allele# original - random alleles are mutated #population[sample(which(str_detect(Genotype, pattern = chromosome)), size = number_mutants, replace = F),# Genotype := str_replace(Genotype, pattern = paste0(chromosome, "_O"), replacement = paste0(chromosome, "_I"))]# current - entire family is mutated on one chromosome, unless the inbreeding allele is on a hemizygous chromosome, then the minimum number of families are mutated to make 5 mutationsif(chromosome !="W"& chromosome !="Y"){ population[Family_ID %in%sample(unique(population[, Family_ID]), size =1), Genotype :=str_replace(Genotype, pattern =paste0(chromosome, "_O"), replacement =paste0(chromosome, "_I"))] } else{# number_mutants is loaded from the parameter table mutations_made <-0# Iterate through unique families as needed Families <-sample(unique(population$Family_ID), Starting_pop_size/N) Family_index <-1while (mutations_made < number_mutants && Family_index <=length(Families)) { current_family <- Families[Family_index]# Identify mutate-able rows from the current family Family_rows <- population[Family_ID == current_family &str_detect(Genotype, pattern = chromosome)]# Calculate how many mutations we can make in this family changes_needed <- number_mutants - mutations_made changes_to_make <-min(nrow(Family_rows), changes_needed)# Select the individuals to mutate rows_to_modify <- Family_rows[1:changes_to_make]# Modify the genotype column in the selected rows population[rows_to_modify, Genotype :=str_replace(Genotype, pattern =paste0(chromosome, "_O"), replacement =paste0(chromosome, "_I")), on ="Individual_ID"]# Update the count of changes made mutations_made <-nrow(population[str_detect(Genotype, "I")])# Move to the next group Family_index <- Family_index +1 } }# Determining the next event# check when the next death occurs: this is the sum of the mortality rates for all individuals in the population next_death <-rexp(n =1, rate =sum(population[, mortality_rate]))# check when the next receptive female-male encounter occurs receptive_females <- population[Sex >0]$Individual_ID # everyone's receptive at the start receptive_males <- population[Sex <1]$Individual_ID# below is the slow way, which I still include to show it's equivalent to the fast way#encounter_possibilities <- # CJ(Female_ID = population[Sex > 0]$Individual_ID,# Male_ID = population[Sex < 1]$Individual_ID)[, encounter_rate := v/number_females]# Find the time the next encounter occurs: plug the sum of the rates into the exponential function. # The population level encounter rate is the product of the rate at which a single male finds a single female, the number of receptive females in the population, and the number of receptive males in the population# next_encounter <-rexp(n =1, rate =length(receptive_females)*length(receptive_males)*v)# Initialize the timer t to the first encounter t <-pmin(next_death, next_encounter)# With the initial population ready to go and the first event found, start the timer and let the simulation run. In short, time progresses as events occur. Events can trigger state changes for the individuals in the population, leading to death, mating and offspring production.while (t <= time_end & keep_going){# death and mating: what type of encounter happens at time tif(next_death < next_encounter){ who_died <- population[sample(.N, 1, prob = mortality_rate)]$Individual_ID # remove individual from population table population <- population[Individual_ID != who_died]# check if mortality event frees up breeding site current_breeders <-sum(population$breeding >0)# If there is an available breeding site, and at least one female to fill it, recruit a new breederif(current_breeders < breeding_sites &&sum(population$Sex >0& population$breeding <1) >0){# assign the new breeders population <- population[sample_vec(which(breeding <1& Sex >0),size =1), # note that all living females have equal prob of becoming a breeder breeding :=1] } } else{ which_encounter <-c(sample_vec(receptive_females, 1), sample_vec(receptive_males, 1)) mates <-NULL# reset this every time as a safeguard# Opposite sex encounters female <- population[Individual_ID %chin% which_encounter[1]] male <- population[Individual_ID %chin% which_encounter[2]]# First determine if a homogametic individual, heterozygous for the I allele, will inbreed on this occasion. heterozygote_inbreeds <-rbinom(1, 1, prob = dominance)# find a brother for the female to encounter, if required brother_encountered <- population[Family_ID %chin% female$Family_ID# assign a receptive sibling - this is the behavioural trait under selection: # should an individual accept an inbred mating early in life occur or not. ][Sex <1& t > refractory_period_end, .SD[sample(.N, 1)]] if(nrow(brother_encountered)>0& female$relative_encountered <1){# do inbreedingif(# female heterozygotes sex_expressed >0&str_detect(female$Genotype, heterozygous_genotype) & heterozygote_inbreeds >0|# female homozygotes sex_expressed >0&str_detect(female$Genotype, homozygous_genotype) |# female hemizygotes sex_expressed >0&str_detect(female$Genotype, hemizygous_genotype) |# male heterozygotes sex_expressed <1&str_detect(brother_encountered$Genotype, heterozygous_genotype) & heterozygote_inbreeds >0|# male homozygotes sex_expressed <1&str_detect(brother_encountered$Genotype, homozygous_genotype) |# male hemizygotes sex_expressed <1&str_detect(brother_encountered$Genotype, hemizygous_genotype)){# if mating occurred, update the population mates <-rbindlist(list(female, brother_encountered)) population[mates,`:=`(matings = matings +1,inbred_mating =fifelse(Sex >0, 1, 0),mated_with =fifelse(Sex >0, brother_encountered$Genotype, NA),refractory_period_end =fifelse(Sex <1, t + (refractory_period * baseline_mean_lifespan), 0)), on = .(Individual_ID)] }# females that had not encountered a relative early in life are coded to have now done so population[female, relative_encountered :=1, on = .(Individual_ID)] }if(female$relative_encountered <1&nrow(brother_encountered) <1){# females that had no receptive brother to encounter are recorded as having had their chance for inbreeding early in life. When the male refractory period != 0, this is unlikely (because all siblings are produced at the same time) but possible. Most commonly, this will occur when a female produces an all-female brood (0.03125 probability) population[female, relative_encountered :=1, on = .(Individual_ID)] }# If the individual has already encountered their sibling, don't swap and let encounter proceed. # If the pair happen to be from the same family, assume individuals can recognise this and depending on genotype, avoid or accept the inbred matingif(# female heterozygotes sex_expressed >0&# female expressed allele female$relative_encountered >0&# already encountered their relative female$Family_ID == male$Family_ID &# they've run into another relative str_detect(female$Genotype, heterozygous_genotype) &# they carry 1 inbreeding allele heterozygote_inbreeds >0|# the allele is expressed# female homozygotes sex_expressed >0&# female expressed allele male$relative_encountered >0&# already encountered relative female$Family_ID == male$Family_ID &# # they've run into another relativestr_detect(female$Genotype, homozygous_genotype) |# female hemizygotes sex_expressed >0& male$relative_encountered >0&# already encountered relative female$Family_ID == male$Family_ID &# # they've run into another relativestr_detect(female$Genotype, hemizygous_genotype) |# male heterozygotes sex_expressed <1&# male expressed allele female$relative_encountered >0&# already encountered relative female$Family_ID == male$Family_ID &# they've encountered another relative str_detect(male$Genotype, heterozygous_genotype) &# they carry 1 inbreeding allele heterozygote_inbreeds >0|# the allele is expressed# male homozygotes sex_expressed <1&# male expressed allele female$relative_encountered >0&# already encountered their relative female$Family_ID == male$Family_ID &# # they've run into another relativestr_detect(male$Genotype, homozygous_genotype) |# male hemizygotes sex_expressed <1&# male expressed allele female$relative_encountered >0&# already encountered their relative female$Family_ID == male$Family_ID &# # they've run into another relativestr_detect(male$Genotype, hemizygous_genotype)){ mates <-rbindlist(list(female, male)) population[mates,`:=`(matings = matings +1,inbred_mating =fifelse(Sex >0, 1, 0),mated_with =fifelse(Sex >0, male$Genotype, NA),refractory_period_end =fifelse(Sex <1, t + (refractory_period * baseline_mean_lifespan), 0)), on = .(Individual_ID)] }# encounters between non-relativesif(female$Family_ID != male$Family_ID & female$relative_encountered >0){ mates <-rbindlist(list(female, male)) population[mates,`:=`(matings = matings +1,mated_with =fifelse(Sex >0, male$Genotype, NA),refractory_period_end =fifelse(Sex <1, t + (refractory_period * baseline_mean_lifespan), 0)), on = .(Individual_ID)] } }# Are there consequences of death and mating: reproduction# check if a female can now produce offspring, either because they're previously mated and have secured a breeding site or because they already held a breeding site and have now mated new_mated_breeder <- population[Sex >0& matings >0& breeding >0& reproduced <1, .(Individual_ID, inbred_mating,Female_genotype = Genotype,Male_genotype = mated_with)]if(nrow(new_mated_breeder) >0) {# add offspring to the population. Each mated female that holds a breeding site produces N offspring offspring <- new_mated_breeder[inheritance_scheme, on = .(Female_genotype = Female_genotype,Male_genotype = Male_genotype), nomatch =NULL, allow.cartesian =TRUE ][, .SD[sample(.N,size = N, prob = zygote_freq, replace = T)] ][, Family_ID := .GRP + Family_ID_counter # assign these offspring to a new family ][, .(Genotype, Sex, inbred_mating, Family_ID) ][, `:=`(mortality_rate =fifelse(inbred_mating >0, 1/(baseline_mean_lifespan + D), # the cost of inbreeding: D <= 01/baseline_mean_lifespan), # outbred offspring mortality ratebirth_time = t,breeding =0,matings =0,reproduced =0,mated_with ="NA",refractory_period_end = t,relative_encountered =0,Individual_ID = .I + Individual_ID_counter,inbred_mating =0)]# bind the offspring table to the existing population table and update which females have reproduced population <-rbindlist(list(population, offspring), use.names =TRUE )[new_mated_breeder, reproduced :=1, on = .(Individual_ID)]# update the Individual_ID counter Individual_ID_counter <-max(population$Individual_ID) Family_ID_counter <-max(population$Family_ID) }# Calculate the frequency of the I allele, quit early if I fixes or goes extinct# calc allele freq if autosomal locus if(chromosome =="A"){ prop_i <- (length(population$Genotype[str_detect(population$Genotype, heterozygous_genotype)]) +2*length(population$Genotype[str_detect(population$Genotype, homozygous_genotype)]))/ (nrow(population)*2) }# calc allele freq if hemizygous locus: W, Y or cytoplasmicif(chromosome =="W"| chromosome =="Y"| chromosome =="C"){ prop_i <- (length(population$Genotype[str_detect(population$Genotype, hemizygous_genotype)])/length(population$Genotype[str_detect(population$Genotype, chromosome)])) }# calc allele freq if diploid in one sex and haploid in the other: X and Zif(chromosome =="X"| chromosome =="Z"){ prop_i <-if(hemizygous_genotype =="X_IY_O"){ (length(population$Genotype[str_detect(population$Genotype, heterozygous_genotype)]) +2*length(population$Genotype[str_detect(population$Genotype, homozygous_genotype)]) +length(population$Genotype[str_detect(population$Genotype, hemizygous_genotype)]))/ (nrow(population[Sex >0])*2+nrow(population[Sex <1]))}else{ (length(population$Genotype[str_detect(population$Genotype, heterozygous_genotype)]) +2*length(population$Genotype[str_detect(population$Genotype, homozygous_genotype)]) +length(population$Genotype[str_detect(population$Genotype, hemizygous_genotype)]))/ (nrow(population[Sex <1])*2+nrow(population[Sex >0]))} }# this is a diagnostic to make sure the model is running well - it can be commented out when running the big simulation prop_i_table <-rbindlist(list(prop_i_table, list(t, prop_i, nrow(population))))print(paste0("Population size = ", nrow(population), ", breeders = ", sum(population$breeding >0), ", time = ", round(t, 3)))if(prop_i >0.5| prop_i <0.0001|nrow(population) <2) keep_going <-FALSE# Move t to next encounter# determining the next event# check when the next death occurs next_death <-rexp(n =1, rate =sum(population[, mortality_rate]))# check when the next receptive female-male encounter occurs receptive_females <- population[Sex >0& matings <1]$Individual_ID receptive_males <- population[Sex <1& refractory_period_end <= t]$Individual_IDif(length(receptive_females)*length(receptive_males) >0){# check the time: this is the sum of the rates, because each male finds each female at the same rate next_encounter <-rexp(n =1, rate =length(receptive_females)*length(receptive_males)*v)# Initialize the timer t to the next encounter }else{next_encounter <-10^3} # make this a number that will never be exceeded by a death time t <- t +pmin(next_death, next_encounter) } finish_time <- t final_pop_size =nrow(population)# cut the i table down to 1 row per 0.1 time increment, join with the relevant parameter space info and save as a csv. prop_i_table %>%mutate(time_group =floor(time /0.1) *0.1) %>%group_by(time_group) %>%slice(1) %>%ungroup() %>%select(-time_group) %>%bind_rows(prop_i_table %>%filter(row_number()==n())) %>%distinct() %>%bind_cols(parameters[row, ] %>%select(-c(heterozygous_genotype, homozygous_genotype, hemizygous_genotype, number_mutants, baseline_mean_lifespan, N, time_end))) %>%write_csv(paste("sim_results/rowID_", parameter_space_ID, chromosome, ".csv", sep =""))}
Define the parameter space
The effect of the relatedness coefficient on the fitness of an inbreeding allele when there is non-zero inbreeding depression is well described. Therefore, we choose to set \(r = 0.5\) for all inbreeding events and instead focus on the effect of varying 1) \(v\), 2) \(\delta\), 3) the male cost to mating, here coded as a refractory period following mating, and 4) the chromosome upon which the the allele is found.
Populations also differ in the number of breeding females they can support, in order to standardise the number of chromosomes carrying the focal locus. For example, when the inbreeding locus is found on an autosome, a population of 200 individuals harbours 400 copies of the inbreeding locus. If the allele is instead found on a cytoplasmic chromosome (and we assume no heteroplasmy), each individual only carries one copy of the locus and a population size of 400 individuals is needed to equalise the population size of chromosomes. \(v\) is adjusted in accordance with the individual population size, so as to ensure that males meet the same number of individuals on average across simulation runs - this is akin to keeping the population density constant between simulation runs.
We run one simulation per row in the parameter space table.
Code
parameters_autosome <- parameters_autosome %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_X <- parameters %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_Z <- parameters_Z %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_Y <- parameters_Y %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_W <- parameters_W %>%slice_sample(prop =1) # shuffle to equalise workload across jobsparameters_C <- parameters_C %>%slice_sample(prop =1) # shuffle to equalise workload across jobs# run the simulation in parallel using all local CPUs on your computer. The code below uses the autosomal parameter space as an example.if(!file.exists("results/complete_results.csv")){# get the number of cores on the computer n.cores <- parallel::detectCores()# make the cluster cluster <- parallel::makeCluster(n.cores)# split parameters across cluster parallel::clusterExport(cluster, paste("parameters_Y")) parallel::parLapply(cluster,1:nrow(parameters_Y), continuous_time_simulation, parameters_Y, offspring_genotypes_Y)}
Load the results
Code
# build a function to load the individual runs and join them into a single tibbleload_results <-function(chromosome){ files <-list.files(path ="sim_results") %>%str_subset(chromosome)paste("sim_results/", files, sep ="") %>%vroom() }if(!file.exists("results/complete_results.csv")){ cytoplasmic_results <-load_results("C") y_sim_results <-load_results("Y") X_sim_results <-load_results("X") Z_sim_results <-load_results("Z") W_sim_results <-load_results("W") A_sim_results <-load_results("A") sim_results <-bind_rows(A_sim_results, y_sim_results, X_sim_results, Z_sim_results, W_sim_results, cytoplasmic_results) sim_results %>%write_csv("results/complete_results.csv")} else{ sim_results <-read_csv("results/complete_results.csv") }
Plot the results
Simulation diagnostics
Code
# load a palette used in all tabstemp <-pnw_palette("Shuksan2",100)
